home *** CD-ROM | disk | FTP | other *** search
/ PC-Blue - MS DOS Public Domain Library / PC-Blue MS-DOS Public Domain Library - NYACC.iso / vol230 / spline.c < prev    next >
Encoding:
C/C++ Source or Header  |  1986-11-30  |  13.0 KB  |  473 lines
  1. /*    spline - interpolate points in a tabulated function
  2.  
  3.     Usage...
  4.         spline [file][options]
  5.         options are:
  6.           -a  [step [start]] automatic abscissas 
  7.           -b                 break interpolation after each label 
  8.           -c                 general curve 
  9.           -n  num            interpolate over num intervals 
  10.           -q                 increase intervals fourfold 
  11.           -t  num            tension in interpolating curve 
  12.           -x  [min [max]]    interpolate from min to max 
  13.           -xl                take logs of x values before interpolating 
  14.           -yl                take logs of y values before interpolating 
  15.  
  16.     Notes...
  17.         This program fits a smooth curve through a given set of points,
  18.         using the splines under tension introduced by Schweikert [J. 
  19.         Math.  and Physics v 45 pp 312-317 (1966)].  They are similar
  20.         to cubic splines, but are less likely to introduce spurious
  21.         inflection points.
  22.  
  23.     History...
  24.         3 Jun 86    -s switch allows specified slopes at ends.
  25.         2 Jun 86    Ver 1.3: Fixed several embarrassing bugs in -b option.
  26.                     -q switch quadruples number of intervals.
  27.         6 Apr 86    -b switch breaks interpolation at labels.
  28.         21 Sep 85    Added -xl and -yl switches for taking logs
  29.         23 Nov 85    Passing lines starting with ';' unchanged, otherwise
  30.                     ignoring them.
  31.  
  32.     Author...
  33.         Copyright (c) 1985, 1986  James R. Van Zandt (jrv@mitre-bedford)
  34.  
  35.         All rights reserved.
  36.         This program may be copied for personal, non-profit use only.
  37.  
  38.         Based on algorithms by A.  K.  Cline [Communications of the ACM
  39.         v 17 n 4 pp 218-223 (Apr 74)].
  40.  
  41. */
  42.  
  43. #include <stdio.h>
  44. #include <math.h>
  45.  
  46. #define VERSION "1.3"
  47. char *RIGHT=        "Copyright (c) 1985, 1986  James R. Van Zandt";
  48.  
  49. #define ENTRIES 200
  50. #define MAXLABELS 50
  51.  
  52. double *x, *y, *temp, *xp, *yp, *path;
  53. double slp1=0.,slpn=0.,sigma=1.;
  54. double abscissa=0.;
  55. double abscissa_step=1.;
  56. double slope1=0.,slopen=0.;    /* slopes supplied by user */
  57. double xmin=0.;
  58. double xmax=0.;
  59. double curv2(), mylog();
  60.  
  61. int xlog=0;        /* nonzero if taking log of x values */
  62. int ylog=0;        /* nonzero if taking log of y values */
  63. int debugging=0;
  64. int breaking=0;        /* nonzero if breaking at labels */
  65. int labels=0;        /* number of labels in input data */
  66. int quadruple=0;    /* nonzero if user asking for 4 times the resolution */
  67. int automatic_abscissas=0;
  68. int slopes=0;        /* nonzero if slopes supplied by user */
  69. int abscissa_arguments=0;
  70. int curve=0;        /* nonzero if general curve permitted */
  71. int x_arguments=0;
  72. int n;                /* number of entries in x, y, yp, and temp */
  73. int index_array[MAXLABELS];    /* indexes into x and y */
  74. int *p_data=index_array;
  75. int total=100;
  76.  
  77. FILE file;
  78.  
  79. char *label_array[MAXLABELS];
  80. char **p_text=label_array;
  81.  
  82. main(argc,argv) int argc; char **argv;
  83. {    int nn,origin;
  84.     read_data(argc,argv);
  85.     if(breaking && labels)
  86.         {origin=0;
  87.         while(labels--)
  88.             {p_data[0] -= origin;
  89.             nn=p_data[0]+1;
  90.             if(quadruple) total=nn*4-3;
  91.             if(nn) curv0(nn,total);
  92.             origin += nn;
  93.             path += nn;
  94.             x += nn;
  95.             y += nn;
  96.             p_data++;
  97.             p_text++;
  98.             }
  99.         }
  100.     else
  101.         {if(quadruple) total=n*4-3;
  102.         curv0(n,total);
  103.         }
  104. }
  105.  
  106. curv0(n,requested)
  107. int n,requested;
  108. {    int i, j, each, stop, seg=0, regressing=0;
  109.     double s,ds,xx,yy;
  110.     for(i=1; i<n; i++) if(path[i]<=path[i-1]) regressing++;
  111.     if (regressing)
  112.         {fprintf(stderr,"ERROR: independent variable not strictly increasing\n");
  113.         return;
  114.         }
  115.     if (curve && (xmin<0. || xmax>1.))
  116.         {fprintf(stderr,"ERROR: independent variable not in range 0 to 1\n");
  117.         return;
  118.         }
  119.     if(curve) {curv1(path,x,xp,n); curv1(path,y,yp,n);}
  120.     else {slp1=slope1; slpn=slopen; curv1(x,y,yp,n);}
  121.     s=path[0];
  122.     seg=0;
  123.     if(requested<n) requested=n;
  124.     if(x_arguments==2)            /* specific range was requested */
  125.         {ds=(xmax-xmin)/requested;
  126.         if(curve)
  127.             {ds*=(path[n-1]-path[0]);
  128.             s=xmin*(path[n-1]-path[0])+path[0];
  129.             }
  130.         else s=xmin;
  131.         for(i=requested+1; i; i--)
  132.             {value(s,&xx,&yy,n);
  133.             printf("%g %g ",xx,yy);
  134.             if(j==p_data[seg]) puts(p_text[seg++]);
  135.             putchar('\n');
  136.             s+=ds;
  137.             }
  138.         }
  139.     else                        /* spline for entire curve was requested */
  140.         {for (i=j=0; j<n-1; j++)
  141.             {stop=requested*(j+1)/(n-1);
  142.             ds=(path[j+1]-path[j])/(stop-i);
  143.             for ( ; i<stop; i++)
  144.                 {value(s,&xx,&yy,n);
  145.                 printf("%g %g ",xx,yy);
  146.                 if(j==p_data[seg]) puts(p_text[seg++]);
  147.                 putchar('\n');
  148.                 s+=ds;
  149.                 }
  150.             }
  151.             xx=x[n-1]; if(xlog) xx=exp(xx);
  152.             yy=y[n-1]; if(ylog) yy=exp(yy);
  153.             printf("%g %g ",xx,yy);
  154.             if(j==p_data[seg]) puts(p_text[seg++]);
  155.             putchar('\n');
  156.         }
  157. }
  158.  
  159. value(s,q,r,n) double s,*q,*r; int n;
  160. {    if(curve)
  161.         {*q=curv2(path,x,xp,s,n);
  162.         *r=curv2(path,y,yp,s,n);
  163.         }
  164.     else
  165.         {*q=s;
  166.         *r=curv2(x,y,yp,s,n);
  167.         }
  168.     if(xlog) *q=exp(*q);
  169.     if(ylog) *r=exp(*r);
  170. }
  171.  
  172. read_data(argc,argv) int argc; char **argv;
  173. {    int i,j,length;
  174.     double xx,yy,d,*pd,sum;
  175.     char *s,*t;
  176. #define BUFSIZE 200
  177.     static char buf[BUFSIZE];
  178.  
  179.     x=path=malloc(ENTRIES*sizeof(double));
  180.     y=malloc(ENTRIES*sizeof(double));
  181.     if(x==0 || y==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  182.     if(argc>1 && streq(argv[1],"?")) help();
  183.     if(argc<=1 || *argv[1]=='-') file=stdin;
  184.     else
  185.         {if(argc>1)
  186.             {file=fopen(argv[1],"r");
  187.             if(file==0) {printf("file %s not found\n",argv[1]); exit();}
  188.             argc--; argv++;
  189.             }
  190.         else help();
  191.         }
  192.     argc--; argv++;
  193.     while(argc>0)
  194.         {i=get_parameter(argc,argv);
  195.         argv=argv+i; argc=argc-i;
  196.         }
  197.     if(sigma<=0.)
  198.         {fprintf(stderr,"ERROR: tension must be positive\n");
  199.         exit(1);
  200.         }
  201.     if(slopes && curve)
  202.         {fprintf(stderr,"ERROR: slopes can\'t be specified for general curve\n");
  203.         exit(1);
  204.         }
  205.     sigma *= -1.;
  206.     if(xlog && !curve)
  207.         {if(x_arguments>1) xmax=mylog(xmax);
  208.         if(x_arguments>=1) xmin=mylog(xmin);
  209.         }
  210.     if(automatic_abscissas && abscissa_arguments<2 && x_arguments>=1)
  211.         abscissa=xmin;
  212.     p_data[0]=-1;
  213.     i=0;
  214.     while(i<ENTRIES)
  215.         {if(fgets(buf,BUFSIZE,file)==0)break;
  216.         t=buf;
  217.         while(*t && isspace(*t)) t++;
  218.         if(*t == 0) continue;        /* skip blank lines */
  219.         buf[strlen(buf)-1]=0;        /* zero out the line feed */
  220.         if(buf[0]==';') {printf("%s\n",buf); continue;}  /* skip comment */
  221.         if(automatic_abscissas)
  222.             {x[i]=abscissa;
  223.             abscissa+=abscissa_step;
  224.             sscanf(buf,"%F",&y[i]);
  225.             if(ylog) y[i]=mylog(y[i]);
  226.             }
  227.         else
  228.             {sscanf(buf,"%F %F",&x[i],&y[i]);
  229.             if(xlog) x[i]=mylog(x[i]);
  230.             if(ylog) y[i]=mylog(y[i]);
  231.             }
  232.         s=buf;                      /* start looking for label */
  233.         while(*s==' ')s++;            /* skip first number */
  234.         while(*s && (*s!=' '))s++;
  235.         if(!automatic_abscissas)    /* skip second number */
  236.             {while(*s==' ')s++;
  237.             while(*s && (*s!=' '))s++;
  238.             }
  239.         while(*s==' ')s++;
  240.         if((length=strlen(s))&&(labels<MAXLABELS))
  241.             {if(*s=='\"')           /* label in quotes */
  242.                 {t=s+1;
  243.                 while(*t && (*t!='\"')) t++;
  244.                 t++;
  245.                 }
  246.             else                    /* label without quotes */
  247.                 {t=s;
  248.                 while(*t && (*t!=' '))t++;
  249.                 }
  250.             *t=0;
  251.             length=t-s;
  252.             p_data[labels]=i;
  253.             p_text[labels]=malloc(length+1);
  254.             if(p_text[labels]) strcpy(p_text[labels++],s);
  255.             }
  256.         i++;
  257.         }
  258.     n=i;
  259.     if(breaking && (!labels || p_data[labels-1]!=n-1))
  260.         {p_data[labels]=i-1;
  261.         if(p_text[labels]=malloc(1)) *p_text[labels++]=0;
  262.         }
  263.     yp=malloc(n*sizeof(double));
  264.     temp=malloc(n*sizeof(double));
  265.     if(temp==0 || yp==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  266.     if(curve)
  267.         {xp=malloc(n*sizeof(double));
  268.         path=malloc(n*sizeof(double));
  269.         if(xp==0|| path==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  270.         path[0]=sum=0.;
  271.         for (i=1; i<n; i++)
  272.             {xx=x[i]-x[i-1];
  273.             yy=y[i]-y[i-1];
  274.             path[i]=(sum+=sqrt(xx*xx + yy*yy));
  275.             }
  276. /*        for(i=0; i<n; i++)
  277.             printf("path[%d]=%g  x[%d]=%g \n",i,path[i],i,x[i]); */
  278.         }
  279. }
  280.  
  281.  
  282. /* get_parameter - process one command line option
  283.         (return # parameters used) */
  284. get_parameter(argc,argv) int argc; char **argv;
  285. {    int i;
  286.     if(streq(*argv,"-a"))
  287.         {i=get_double(argc,argv,2,&abscissa_step,&abscissa,&abscissa);
  288.         abscissa_arguments=i-1;
  289.         automatic_abscissas=1;
  290.         return i;
  291.         }
  292.     else if(streq(*argv,"-b")) {breaking=1; return 1;}
  293.     else if(streq(*argv,"-c")) {curve=1; return 1;}
  294.     else if(streq(*argv,"-d")) {debugging=1; return 1;}
  295.     else if(streq(*argv,"-n"))
  296.         {if((argc>1) && numeric(argv[1])) total=atoi(argv[1]);
  297.         return 2;
  298.         }
  299.     else if(streq(*argv,"-q")) {quadruple=1; return 1;}
  300.     else if(streq(*argv,"-t"))
  301.         {return(get_double(argc,argv,1,&sigma,&abscissa,&abscissa));
  302.         }
  303.     else if(streq(*argv,"-s"))
  304.         {slopes++;
  305.         return(get_double(argc,argv,2,&slope1,&slopen,&slopen));
  306.         }
  307.     else if(streq(*argv,"-x"))
  308.         {i=get_double(argc,argv,2,&xmin,&xmax,&xmax);
  309.         x_arguments=i-1;
  310.         return i;
  311.         }
  312.     else if(streq(*argv,"-xl")) {xlog++; return 1;}
  313.     else if(streq(*argv,"-yl")) {ylog++; return 1;}
  314.     else gripe(argv);
  315. }
  316.  
  317. get_double(argc,argv,permitted,a,b,c)
  318. int argc,permitted; char **argv; double *a,*b,*c;
  319. {    int i=1;
  320.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *a=atof(argv[i++]);
  321.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *b=atof(argv[i++]);
  322.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *c=atof(argv[i++]);
  323.     return i;
  324. }
  325.  
  326. int streq(a,b) char *a,*b;
  327. {    while(*a)
  328.         {if(*a!=*b)return 0;
  329.         a++; b++;
  330.         }
  331.     return 1;
  332. }
  333.  
  334. gripe_arg(s) char *s;
  335. {    fprintf(stderr,"argument missing for switch %s",s);
  336.     help();
  337. }
  338.  
  339. gripe(argv) char **argv;
  340. {    fprintf(stderr,*argv); fprintf(stderr," isn\'t a legal argument \n\n");
  341.     help();
  342. }
  343.  
  344. help()
  345. {    fprintf(stderr,"spline   version %s",VERSION);
  346.     fprintf(stderr,"\nusage: spline [file][options]\n");
  347.     fprintf(stderr,"options are:\n");
  348.     fprintf(stderr,"  -a  [step [start]] automatic abscissas \n");
  349.     fprintf(stderr,"  -b                 break interpolation after each label \n");
  350.     fprintf(stderr,"  -c                 general curve \n");
  351.     fprintf(stderr,"  -n  num            interpolate over num intervals \n");
  352.     fprintf(stderr,"  -q                 increase intervals fourfold \n");
  353.     fprintf(stderr,"  -s  [num [num]]    specify slopes at beginning and end \n");
  354.     fprintf(stderr,"  -t  num            tension in interpolating curve \n");
  355.     fprintf(stderr,"  -x  [min [max]]    interpolate from min to max \n");
  356.     fprintf(stderr,"  -xl                take logs of x values before interpolating \n");
  357.     fprintf(stderr,"  -yl                take logs of y values before interpolating \n");
  358.     exit();
  359. }
  360.  
  361. numeric(s) char *s;
  362. {    char c;
  363.     while(c=*s++)
  364.         {if((c<='9' && c>='0') || c=='+' || c=='-' || c=='.') continue;
  365.         return 0;
  366.         }
  367.     return 1;
  368. }
  369.  
  370. curv1(x,y,yp,n) double *x,*y,*yp; int n;
  371. {    int i;
  372.     double c1,c2,c3,deln,delnm1,delnn,dels,delx1,delx2,delx12;
  373.     double diag1,diag2,diagin,dx1,dx2,exps;
  374.     double sigmap,sinhs,sinhin,slpp1,slppn,spdiag;
  375.  
  376.     delx1=x[1] - x[0];
  377.     dx1=(y[1] - y[0])/delx1;
  378.     if(slopes) {slpp1=slp1; slppn=slpn;}
  379.     else
  380.         {if(n!=2)
  381.             {delx2= x[2] - x[1];
  382.             delx12= x[2] - x[0];
  383.             c1= -(delx12 + delx1)/delx12/delx1;
  384.             c2= delx12/delx1/delx2;
  385.             c3= -delx1/delx12/delx2;
  386.             slpp1= c1*y[0] + c2*y[1] + c3*y[2];
  387.             deln= x[n-1] - x[n-2];
  388.             delnm1= x[n-2] - x[n-3];
  389.             delnn= x[n-1] - x[n-3];
  390.             c1= (delnn + deln)/delnn/deln;
  391.             c2= -delnn/deln/delnm1;
  392.             c3= deln/delnn/delnm1;
  393.             slppn= c3*y[n-3] + c2*y[n-2] + c1*y[n-1];
  394.             }
  395.         else yp[0]=yp[1]=0.;
  396.         }
  397.                         /* denormalize tension factor */
  398.     sigmap=fabs(sigma)*(n-1)/(x[n-1]-x[0]);
  399.                         /* set up right hand side and tridiagonal system for
  400.                            yp and perform forward elimination                */
  401.     dels=sigmap*delx1;
  402.     exps=exp(dels);
  403.     sinhs=.5*(exps-1./exps);
  404.     sinhin=1./(delx1*sinhs);
  405.     diag1=sinhin*(dels*.5*(exps+1./exps)-sinhs);
  406.     diagin=1./diag1;
  407.     yp[0]=diagin*(dx1-slpp1);
  408.     spdiag=sinhin*(sinhs-dels);
  409.     temp[0]=diagin*spdiag;
  410.     if(n!=2)
  411.         {for(i=1; i<=n-2; i++)
  412.             {delx2=x[i+1]-x[i];
  413.             dx2=(y[i+1]-y[i])/delx2;
  414.             dels=sigmap*delx2;
  415.             exps=exp(dels);
  416.             sinhs=.5*(exps-1./exps);
  417.             sinhin=1./(delx2*sinhs);
  418.             diag2=sinhin*(dels*(.5*(exps+1./exps))-sinhs);
  419.             diagin=1./(diag1+diag2-spdiag*temp[i-1]);
  420.             yp[i]=diagin*(dx2-dx1-spdiag*yp[i-1]);
  421.             spdiag=sinhin*(sinhs-dels);
  422.             temp[i]=diagin*spdiag;
  423.             dx1=dx2;
  424.             diag1=diag2;
  425.             }
  426.         }
  427.     diagin=1./(diag1-spdiag*temp[n-2]);
  428.     yp[n-1]=diagin*(slppn-dx2-spdiag*yp[n-2]);
  429.                     /* perform back substitution */
  430.     for (i=n-2; i>=0; i--) yp[i] -= temp[i]*yp[i+1];
  431. }
  432.  
  433.  
  434. double curv2(x,y,yp,t,n) double *x,*y,*yp,t; int n;
  435. {    int i,j;
  436.     static int i1=1;
  437.     double del1,del2,dels,exps,exps1,s,sigmap,sinhd1,sinhd2,sinhs;
  438.  
  439.     s=x[n-1]-x[0];
  440.     sigmap=fabs(sigma)*(n-1)/s;
  441. #ifdef WORK
  442.     for (j=2; j; j--)        /* want: x[i-1] <= t < x[i], 0 < i <= n */
  443.         {for (i=i1; i<n; i++)
  444.             {if(x[i]>t) break;
  445.             }
  446.         if(i==n) i=n-1;
  447.         if(x[i-1]<=t || t<=x[0]) break;
  448.         i1=1;
  449.         }
  450. #endif
  451.     i=i1;
  452.     if(i>n) i=1;
  453.     while(i<n && t>=x[i]) i++;
  454.     while(i>1 && x[i-1]>t) i--;
  455.     i1=i;
  456.     del1=t-x[i-1];
  457.     del2=x[i]-t;
  458.     dels=x[i] - x[i-1];
  459.     exps1=exp(sigmap*del1); sinhd1=.5*(exps1-1./exps1);
  460.     exps= exp(sigmap*del2); sinhd2=.5*(exps-1./exps);
  461.     exps=exps1*exps;        sinhs=.5*(exps-1./exps);
  462.     return ((yp[i]*sinhd1 + yp[i-1]*sinhd2)/sinhs +
  463.             ((y[i] - yp[i])*del1 + (y[i-1] - yp[i-1])*del2)/dels);
  464. }
  465.  
  466. double mylog(x) double x;
  467. {    if(x>0.) return log(x);
  468.     fprintf(stderr,"can%'t take log of nonpositive number");
  469.     exit(1);
  470.     return 0.;
  471. }
  472.  
  473.